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Abstract 



O 

D . Identification of time-varying linear systems, which introduce both time-shifts (delays) and frequency- 



shifts (Doppler-shifts), is a central task in many engineering applications. This paper studies the problem 
of identification of underspread linear systems (ULSs), whose responses lie within a unit-area region in the 
delay-Doppler space, by probing them with a known input signal. It is shown that sufficiently-underspread 
parametric linear systems, described by a finite set of delays and Doppler-shifts, are identifiable from a 
single observation as long as the time-bandwidth product of the input signal is proportional to the square 
of the total number of delay-Doppler pairs in the system. In addition, an algorithm is developed that 
enables identification of parametric ULSs from an input train of pulses in polynomial time by exploiting 

t-h ; 

, recent results on sub-Nyquist sampling for time delay estimation and classical results on recovery of 

00 

frequencies from a sum of complex exponentials. Finally, application of these results to super-resolution 
target detection using radar is discussed. Specifically, it is shown that the proposed procedure allows 

o ; 

to distinguish between multiple targets with very close proximity in the delay-Doppler space, resulting 



in a resolution that substantially exceeds that of standard matched-filtering based techniques without 
introducing leakage effects inherent in recently proposed compressed sensing-based radar methods. 
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Fig. 1. Schematic representation of identification of a time- varying linear system H by probing it with a known input signal. 
Characterization of an identification scheme involves specification of the input probe, x(t), and the accompanying sampling and 
recovery stages. 



I. Introduction 

Physical systems arising in a number of application areas can often be described as linear and time 
varying HI, [H. Identification of such systems may help improve overall performance, e.g., the bit-error 
rate in communications (H, or constitute an integral part of the overall system operation, e.g., target 
detection using radar or active sonar O. 

Mathematically, identification of a given time-varying linear system T~L involves probing it with a known 
input signal x(t) and identifying 7~L by analyzing the single system output l~L{x(t)) 0, as illustrated in 
Fig. [T] Unlike time-invariant linear systems, however, a single observation of a time-varying linear system 
does not lead to a unique solution unless additional constraints on the system response are imposed. This 
is due to the fact that such systems introduce both time-shifts {delays) and frequency-shifts (Doppler- 
shifts) to the input signal. It is now a well-established fact in the literature that a time-varying linear 
system T~L can only be identified from a single observation if l~L(5{t)) is known to lie within a region 1Z 
of the delay-Doppler space such that area(T^) < 1 O-O. Identifiable time-varying linear systems are 
termed underspread, as opposed to nonidentifiable overspread linear systems, which satisfy area(7£) > 1 

0, ml] 

In this paper, we study the problem of identification of underspread linear systems (ULSs) whose 
responses can be described by a finite set of delays and Doppler-shifts. That is, 

K 

n(x(t)) = ^a k x(t-T k )e^ t (1) 

k=l 

where (tj,, v k ) denotes a delay-Doppler pair and a k G C is the complex attenuation factor associated with 
(Tfci Vk). Unlike most of the existing work in the literature, however, our goal in this paper is to explicitly 

'it is still an open research question as to whether critically-spread linear systems, which correspond to are&(7V) = 1, are 
identifiable or nonidentifiable (6); see (5) for a partial answer to this question for the case when 1Z is a rectangular region. 
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characterize conditions on the bandwidth and temporal support of the input signal that ensure identification 
of such ULSs from single observations. The importance of this goal can be best put into perspective 
by realizing that ULSs of the form £[) tend to arise frequently in many applications. Consider, for 
example, a single-antenna transmitter communicating wirelessly with a single-antenna receiver in a mobile 
environment. Over a small-enough time interval, the channel between the transmitter and receiver is of the 
form (Q]) with each triplet (ife, z^, corresponding to a distinct physical path between the transmitter 
and receiver Q. Identification of % can enable one to establish a relatively error-free communication link 
between the transmitter and receiver. But wireless systems also need to identify channels using signals 
that have as small time-bandwidth product as possible so that they can allocate the rest of the temporal 
degrees of freedom to communicating data fTl, l8i 

Similarly, in the case of target detection using radar or active sonar, the (noiseless, clutter-free) received 
signal is of the form (Q]) with each triplet (jk, u^, a*,) corresponding to an echo of the transmitted signal 
from a distinct target in the delay-Doppler space 0. Identification of % in this case enables one to 
accurately obtain the radial position and velocity of the targets. Radar systems also strive to operate with 
signals (waveforms) that have as small temporal support and bandwidth as possible. This is because the 
temporal support of the radar waveform is directly tied to the time it takes to identify all the targets while 
the bandwidth of the waveform — among other technical considerations — is tied to the sampling rate of 
the radar receiver O. 

Given the ubiquity of time-varying linear systems in engineering applications, there exists considerable 
amount of existing literature that studies identification of such systems in an abstract setting. Kailath was 
the first to recognize that the identifiability of a time-varying linear system % from a single observation 
is directly tied to the area of the region 1Z that contains T-L{5{t)) 01. Kailath's seminal work in J4] laid 
the foundations for the future works of Bello (5), Kozek and Pfander (3), and Pfander and Walnut |fj), 
which establish the nonidentifiability of overspread linear systems and provide constructive proofs for 
the identifiability of arbitrary ULSs. However, none of these results shed any light on the bandwidth and 
temporal support of the input signal needed to ensure identification of ULSs of the form £[]). On the 
contrary, the constructive proofs provided in (3l-@ require use of input signals with infinite bandwidth 
and temporal support. 

In contrast, to the best of our knowledge, this is the first paper to develop a theory for identification 
of ULSs of the form (Q]), henceforth referred to as parametric ULSs, that parallels that of O-O for 
identification of arbitrary ULSs. One of the main contributions of this paper is that we establish using a 
constructive proof that sufficiently-underspread parametric linear systems are identifiable as long as the 
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time-bandwidth product of the input signal is proportional to the square of the total number of delay- 
Doppler pairs in the system. Equally importantly, as part of our constructive proof, we concretely specify 
the nature of the input signal (a finite train of pulses) and the structure of a corresponding polynomial- 
time (in the number of delay-Doppler pairs) recovery procedure that enable identification of parametric 
ULSs. These ideas are also immediately applicable to super-resolution target detection using radar and 
we show in the latter part of the paper that our approach indeed results in a resolution that substantially 
exceeds that of standard matched-filtering based techniques without introducing leakage effects inherent 
in recently proposed compressed sensing-based radar methods 0. 

The key developments in the paper leverage recent results on sub-Nyquist sampling for time-delay 
estimation |[T0l and classical results on direction-of-arrival (DOA) estimation lflTl - lfl4l . Unlike the 
traditional DOA estimation literature, however, we do not assume that the system output is observed 
by an array of antennas. Additionally, in contrast to ifTOl . our goal here is not a reduction in the sampling 
rate; rather, we are interested in characterizing the minimum temporal degrees of freedom of the input 
signal needed to ensure identification of a parametric ULS %. The connection to sub-Nyquist sampling 
can be understood by noting that the sub-Nyquist sampling results of IfTOl enable recovery of the delays 
associated with % using a small-bandwidth input signal. Further, the "train-of-pulses" nature of the input 
signal combined with results on recovery of frequencies from a sum of complex exponentials Ifl4l allow 
recovery of the Doppler-shifts and attenuation factors using an input signal of small temporal support. 

Several works in the past have considered identification of specialized versions of parametric ULSs. 
Specifically, (9], |[T5l - |[T8l treat parametric ULSs whose delays and Doppler-shifts lie on a quantized 
grid in the delay-Doppler space. On the other hand, |[T9l considers the case in which there are no more 
than two Doppler-shifts associated with the same delay. The proposed recovery procedures in |[T9l also 
have exponential complexity, since they require exhaustive searches in a A'-dimensional space, and stable 
initializations of these procedures stipulate that the system output be observed by an M-element antenna 
array with M <L K. 

While the insights of 0, lfT5l - |[T8l can be extended to arbitrary parametric ULSs by taking infmitesimally- 
fme quantization of the delay-Doppler space, this will require input signals with infinite bandwidth and 
temporal support. In contrast, our ability to avoid quantization of the delay-Doppler space is due to 
the fact that we treat the system-identification problem directly in the analog domain. This follows the 
philosophy in much of the recent work in analog compressed sensing, termed Xampling, which provides 
a framework for incorporating and exploiting structure in analog signals without the need for quantization 
l|20l - ||25l . This is in particular the key enabling factor that helps us avoid the catastrophic implications 
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of the leakage effects in the context of radar target detection. 

Before concluding this discussion, we note that responses of arbitrary ULSs can always be represented 
as £[]) under the limit K — > oo. Therefore, the main result of this paper can also be construed as an 
alternate constructive proof of the identifiability of sufficiently-underspread linear systems. Nevertheless, 
just like J3l-||6l, this interpretation of the presented results again seem to suggest that identification of 
arbitrary ULSs requires use of input signals with infinite bandwidth and temporal support. 

The rest of this paper is organized as follows. In Section [TTJ we formalize the problem of identification 
of parametric ULSs along with the accompanying assumptions. In Section [Till we propose a polynomial- 
time recovery procedure used for the identification of parametric ULSs, while Section [TV] specifies the 
conditions on the input signal needed to guarantee unique identification using the proposed procedure. 
We compare the results of this paper to some of the related literature on identification of parametric 
ULSs in Section [V] and discuss an application of our results to super-resolution target detection using 
radar in Section [VTJ Finally, we present results of some numerical experiments in Section IVIII 

We make use of the following notational convention throughout this paper. Vectors and matrices are 
denoted by bold-faced lowercase and bold-faced uppercase letters, respectively. The nth element of a 
vector a is written as a n , and the (i, j)th element of a matrix A is denoted by A^. Superscripts (•)*, 
(•) T and (•) represent conjugation, transposition, and conjugate transposition, respectively. In addition, 
the Fourier transform of a continuous-time signal x (t) G L2OC) is defined by X (to) = x (t) e~^ t dt, 
while (x (t) ,y(t)) = f^°x(t)y*(t)dt denotes the inner product between two continuous-time signals 
in L2OC). Similarly, the discrete-time Fourier transform of a sequence a[n] G ^(C) is defined by 
nez a M e ^ nT anc * i s periodic in to with period 2tt/T. Finally, we use A^ to write the 
Moore-Penrose pseudoinverse of a matrix A. 

II. Problem Formulation and Main Results 

In this section, we formalize the problem of identification of a parametric ULS % whose response is 
described by a total of K arbitrary delay-Doppler-shifts of the input signal. The task of identification of 
% essentially requires specifying two distinct but highly intertwined steps. First, we need to specify the 
conditions on the bandwidth and temporal support of the input signal x(t) that ensure identification of % 
from a single observation. Second, we need to provide a polynomial-time recovery procedure that takes 
as input %{x{t)) and provides an estimate % of the system response by exploiting the properties of x(t) 
specified in the first step. We begin by detailing our system model and the accompanying assumptions. 

In (Q]), some of the delays, r^, might be repeated. Expressing % in terms of K T < K distinct delays 
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in this case leads to 

U{x(t)) = Y^Y. a v x{ t ~ ^ w (2) 
i=i j=i 

where z/y denotes the jth Doppler-shift associated with the zth distinct delay Tj, a>ij £ C denotes the 
attenuation factor associated with the delay-Doppler pair (rj,^), and K = Yli=i^v,i- We choose to 
express 1-L(x(t)) as in (f2]) so as to facilitate the forthcoming analysis. Throughout the rest of the paper, 
we use t = {tj, i = 1, . . . ,K T } to denote the set of K T distinct delays associated with T~L. The first 
main assumption that we make concerns the footprint of H in the delay-Doppler space: 
[Al] The response l~L{5{t)) of 7i lies within a rectangular region of the delay-Doppler space; in other 
words, (n,Vij) G [0,r max ] x [-u max /2,u max /2]. This is indeed the case in many engineering 
applications (see, e.g., CQ, Gl)- The parameters T max and v max are termed in the parlance of linear 
systems as the delay spread and the Doppler spread of the system, respectively. 
Next, we use T and W to denote the temporal support and the two-sided bandwidth of the known 
input signal x(t) used to probe T~L, respectively. We propose using input signals that correspond to a finite 
train of pulses: 



N-l 
n=0 



(t) = J2 x n9(t-nT),0<t<T (3) 



where g(t) is a prototype pulse of bandwidth W that is (essentially) temporally supported on [0, T] and 
is assumed to have unit energy (f \g(t)\ 2 dt = 1), and {x n G C} is an iV-length probing sequence. 
The parameter N is proportional to the time-bandwidth product of x(t), which roughly defines the 
number of temporal degrees of freedom available for estimating % [ Sj : N = T/T oc TW^ The final 
two assumptions that we make concern the relationship between the delay spread T max and the Doppler 
spread v max of H, and the temporal support T and bandwidth W of g(t): 

[A2] The delay spread of 7~L is strictly smaller than the temporal support of g(t): r max < T, and 
[A3] The Doppler spread of T~L is much smaller than the bandwidth of g(t): v max <^ W. 

Note that, since W oc 1/T, [A3] equivalently imposes that v m axT <C 1- This assumption states that 
the temporal scale of variations in T-L is large relative to the temporal scale of variations in x(t). It is 
worth pointing out that linear systems that are sufficiently underspread in the sense that T max u max <C 1 
can always be made to satisfy [A2] and [A3] for any given budget of the time-bandwidth product. 

2 Recall that the temporal support and the bandwidth of an arbitrary pulse g(t) are related to each other as W oc 1/T. 
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Remark 1. In order to elaborate on the validity of [A2] and [A3], note that there exist many communication 
applications where underlying linear systems tend to be highly underspread HI § 14.2]. Similarly, [A2] 
and [A3] in the context of radar target detection simply mean that the targets are not too far away from 
the radar and that their velocities are not too high. Consider, for example, an L-band radar (operating 
frequency of 1.3 GHz) that transmits a pulse g(t) of bandwidth W = 10 MHz after every T = 50 //s. 
Then both [A2] and [A3] are satisfied when the distance between the radar and any target is at most 7.5 
km [r max k, 50 ps) and the radial velocity of any target is at most 185 km/h iy max « 445 Hz) ||2). 

The following theorem summarizes our key result concerning identification of parametric ULSs. 

Theorem 1 (Identification of Parametric Underspread Linear Systems). Suppose that % is a parametric 
ULS that is completely described by a total of K = ^^Z\^vi triplets (ri,Vij,a.ij). Then, irrespective 
of the distribution o/{(Tj,i/y)} within the delay— Doppler space, H can be identified in polynomial-time 
from a single observation 7i(x(t)) as long as it satisfies [A1]—[A3], the probing sequence {x n } remains 
bounded away from zero in the sense that \x n \ > Ofor every n = 0, . . . , N — 1, and the time-bandwidth 
product of the known input signal x(t) satisfies the condition 

TW > 8irK T K Utmax (4) 

where K Vjmax = maxj K v i is the maximum number of Doppler-shifts associated with any one of the 
distinct delays. In addition, the time-bandwidth product of x(t) is guaranteed to satisfy © as long as 
TW > 2it{K + l) 2 . 

The rest of this paper is devoted to providing a proof of Theorem Q] In terms of a general roadmap for the 
proof, we first exploit the sub-Nyquist sampling results of iflOl to argue that x(t) with small bandwidth 
suffices to recover the delays associated with %. We then exploit the "train-of-pulses" structure of x(t) and 
classical results on recovery of frequencies from a sum of complex exponentials 11141 to argue that x(t) 
with small temporal support suffices to recover the Doppler-shifts and attenuation factors associated with 
Ti. The statement of Theorem Q] will then follow by a simple combination of the two claims concerning 
the bandwidth and temporal support of x(t). We will make use of © and ([3]) in the following to describe: 
[1] The polynomial-time recovery procedure used for the identification of % (cf. Section UTTT ). and 
[2] The accompanying conditions on x(t) needed to guarantee identification of % (cf. Section ITVT). 
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Fig. 2. Schematic representation of the polynomial-time recovery procedure for identification of parametric underspread linear 
systems from single observations. 



III. Polynomial-Time Identification of ULSs 

In this section, we characterize the polynomial-time recovery procedure used for identification of 
ULSs of the form In order to facilitate understanding of the proposed algorithm, shown in Fig. |2j 
we conceptually partition the method into two stages: sampling and recovery. The rest of this section is 
devoted to describing these two steps in detail. Before proceeding further, however, it is instructive to 
first make use of © and © and rewrite the output of H as 

K T K Uyi N-l 

n(x(t)) = E E E Oijx^^git -n- nT) 

1 = 1 j = l 71=0 

^ E E E CHjXne? 2 ^ nT g(t -n- nT) 

i=l j=l n=0 
K T N-l 

= EE a *M^-Ti-nT) (5) 

i=l n =0 

where (a) follows from the assumption v max T <C 1, which implies that e i 27r ^j* e j 2v ^i nT for all 
t G [(n — 1)T, nT), and the sequences {a,;[n]}, i = 1, . . . , K T , are defined as 

<n[n] = ^2 «■,■<■>■< " '• n = 0, . . . , N - 1. (6) 

3=1 

A. The Sampling Stage 

We leverage the ideas of iflOl on time-delay estimation from sub-Nyquist samples to describe the 
sampling stage of our recovery procedure. While the primary objective in IflOl is time-delay estimation 
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from low-rate samples, the development here is carried out with an eye towards identification of parametric 
ULSs regardless of the distribution of system parameters within the delay-Doppler space — the so- 
called super-resolution identification. In ifTOl , a general multi-channel sub-Nyquist sampling scheme was 
introduced for the purpose of recovering a set of unknown delays from signals of the form ©. Here, 
we focus on one special case of that scheme, which consists of a low-pass filter (LPF) followed by a 
uniform sampler. This architecture may be preferable from an implementation viewpoint since it requires 
only one sampling channel, thereby simplifying analog front-end of the sampling hardware. The LPF, 
besides being required by the sampling stage, also serves as the front-end of the system-identification 
hardware and rejects noise and interference outside the working spectral band. 

Our sampling stage first passes the system output y(t) = %{x{t)) through a LPF whose impulse 
response is given by s*(—t) and then uniformly samples the LPF output at times {t = nT/p}. We 
assume that the frequency response, S*(oj), of the LPF is contained in the spectral band F, defined as 



F 



7T 7T 
-P, ~p 



(7) 



and is zero for frequencies uj ^ F. Here, the parameter p is assumed to be even and satisfies the condition 
p > 2K T ; exact requirements on p to ensure identification of % will be given in Section [TV] In order to 
relate the sampled output of the LPF with the multi-channel sampling formulation of ifTOll . we define p 
sampling (sub)sequences as 

c e [n] = (y(t),s(t-nT-(£-l)T/p)), 1=1,..., p. (8) 

These sequences correspond to periodically splitting the samples at the output of the LPF, which is 
generated at a rate of p/T, into p slower sequences at a rate of 1/T each using a serial-to-parallel 
converter; see Fig. [2] for a schematic representation of this splitting. 

Next, we define the vector c (e Ja;T ) as the p-length vector whose £th element is Cg (e JwT ) , which 
denotes the discrete-time Fourier transform (DTFT) of C([n}. In a similar fashion, we define a (e ju)T ) 
as the ir r -length vector whose ith element is given by Ai (e JajT ), the DTFT of a,i[n]. It can be shown 
following the developments carried out in IfTOl that these two vectors are related to each other by 

c (e^ T ) = W (e^ T ) N (r) D (e^ T , r) a (e^ T ) . (9) 
Here, W (e JwT ) is a p x p matrix with (I, m)th element 

W em (e^ T ) = e ^-DT/ P * <?* L + G r + ^ ^p-xx (1Q) 
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where m' = m — p/2 — 1, N (r) isapx K T Vandermonde matrix with (m, i)th element 

N mi (T) = e-^ TO ' T *, (11) 

andD(e^ T ,r) is a iif T x K T diagonal matrix whose ith diagonal element is given by e ]U)T> . Assuming 
for the time being that W (e JwT ) is a stably-invertible matrix, we define the modified measurement vector 
d (e jujT ) = W" 1 (e jujT ) c (e jujT ). Denoting 

b(e^ T ) =D(e^ T ,r)a( e ^ T ), (12) 

we see from ( fTOl ) that 

d (e jujT ) = N (t) b (e*" r ) . (13) 

Since N (r) is not a function of oj, (TT3T > can be expressed in the discrete-time domain using the linearity 
of the DTFT as 

d[n] = N(r)b[n], neZ, (14) 

Here, the elements of the vectors d[n] and h[n] are discrete-time sequences that are given by the inverse 
DTFT of the elements of d (e JwT ) and b (e^ 7 ) , respectively. 

The key insight to be drawn here is that ([T3l . and its time-domain equivalent ([141 . can be viewed as an 
infinite ensemble of modified measurement vectors in which each element corresponds to a distinct matrix 
N (r) that, in turn, depends on the set of (distinct) delays r. Linear measurement models of the form 
(fT4l — in which the measurement matrix is completely determined by a set of (unknown) parameters — have 
been studied extensively in a number of research areas such as system identification |[26l and direction- 
of-arrival and spectrum estimation lfl4l . |[27l . One specific class of methods that has proven to be quite 
useful in these areas in efficiently recovering the parameters that describe the measurement matrix are 
the so-called subspace methods 1271 . Consequently, our approach in the recovery stage will be to first use 
subspace methods in order to recover the set r from d[n]. Afterwards, since ]NP (r) N (r) = I because 
of the assumption that p > 2K T , we will recover a (e- 7 ^) from d[n] using linear filtering operations as 
follows [cf. (H3, ([13] 

a (e^ T ) = D" 1 (e^ r , r) N+ (r) d (e^ T ) . (15) 

Finally, the Doppler-shifts and attenuation factors associated with H are determined from the vector 
a (e.i wT ) by an additional use of the subspace methods. 
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Before proceeding to the recovery stage, we justify the assumption that W (e Jwr ) can be stably 
inverted. To this end, observe from (fTOb that W (e^ 7 ) can be decomposed as 

W (e^ T ) = * (e jujT ) F H * (e^ T ) , (16) 

where * (e jujT ) is a p x p diagonal matrix with Ah diagonal element 

= y/ft-l^efvV-W*, (17) 

F is a p-point discrete Fourier transform (DFT) matrix with (£, m)th element equal to 

F em = ±e->T^^\ (18) 

and * (e jujT ) is a p x p diagonal matrix whose mth diagonal element is given by 

(e^ T ) = ^S* L + ^(m- p/2 - 1)) G L + ^(m -p/2 - 1)) . (19) 

It can now be easily seen from the decomposition in (fT6l) that, in order for W (e Ja;T ) to be stably 
invertible, each of the above three matrices has to be stably invertible. By construction, both (e J ^ T ) 
and are stably invertible. The invertibility requirement on the diagonal matrix \I/ (e ja;T ) leads to the 
following conditions on the pulse g(t) and the kernel s*(—t) of the LPF. 

Condition 1. In order for \I/ (e jajT ) to be stably invertible, the continuous-time Fourier transform of 
g (t) has to satisfy 

a<\G(u;)\<b a.e.ueJ 7 (20) 
for some positive constants a > and b < oo. 

Condition 2. In order for \l/ (e JwT ) to be stably invertible, the continuous-time Fourier transform of the 
LPF s* (— t) has to satisfy 

c < \S < d a.e.coeJ 7 (21) 
for some positive constants c > and d < oo. 

Condition Q] requires that the bandwidth W of the prototype pulse g(t) has to satisfy 

W > ^p. (22) 

In Section [TV] we will derive additional conditions on the parameter p and combine them with (l22l to 
obtain equivalent requirements on the time-bandwidth product of the input signal x(t) that will ensure 
invertibility of the matrix W (e^ 7 ). 
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We conclude discussion of the sampling stage by pointing out that the decomposition in (fT6l ) also 
provides an efficient way to implement the digital-correction filter bank W" 1 (e- 7 '^). This is because 
( fT6l ) implies that 

W -l (gjVT) = ^-1 (gje-T) F$ -l _ (23) 

Therefore the implementation of W _1 (e JwT ) can be done in three stages, where each stage corresponds 
to one of the three matrices in (|23T ). Specifically, define the set of digital filters {<^[n]} and {^[n]} as 

<j> t [n] = IDTFT { (e^ T ) } [n] , 1 < t < p (24) 

and 

i> t [n] = IDTFT { (e^ T ) } [n] , 1 < £ < p, (25) 

where IDTFT denotes the inverse DTFT operation. The first correction stage involves filtering the 
sequences {c^[n]} using the set of filters {^[n]}. Next, multiplication with the DFT matrix F can 
be efficiently implemented by applying the fast Fourier transform (FFT) to the outputs of the filters in 
the first stage. Finally, the third correction stage involves filtering the resulting sequences using the set 
of filters {^[n]} to get the desired sequences {^[n]}. This last correction stage can be interpreted as an 
equalization step that compensates for the non-flatness of the frequency responses of the prototype pulse 
and the kernel of the LPF A detailed schematic representation of the sampling stage, which is based on 
the preceding interpretation of W" 1 (e Ja;T ), is provided in Fig. [2] 

B. The Recovery Stage 

We conclude this section by describing in detail the recovery stage, which — as noted earlier — consists 
of two steps. In the first step, we rely on subspace methods to recover the delays r from d[n] [cf. (fT4l)l. 
In the second step, we make use of the recovered delays to obtain the Doppler-shifts and attenuation 
factors associated with each of the delays. 

1) Recovery of the Delays: In order to recover r from d[n], we rely on the approach advocated in 
ifTOll and make use of the well-known ESPRIT algorithm |[28l together with an additional smoothing 
stage ||29l . The exact algorithm is given in Table Jl we refer the reader to IfTOll . |[28l for details. 

2 ) Recovery of the Doppler-Shifts and Attenuation Factors: Once the unknown delays are found, we 
can recover the vectors a[n] through the frequency relation (H~5b . Next, define for each delay n, the set 
of corresponding Doppler-shifts 

f i = {vij, j = l,..., K Vti } (26) 
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TABLE I 

Delay Recovery Algorithm 



(i) Construct the matrix 

— 1 M 

where d m is the M = p/2 length subvector which is given by 

d m [n] = ^ d m [n] d m +i [n] ... d m+A/ [n] j • 

(ii) Recover K T as the rank of Rdd- 

(iii) Perform a singular value decomposition (SVD) of Udd and construct the matrix E s consisting of the K T singular vectors 
corresponding to the K T nonzero singular values of Rdd as its columns. 

(iv) Compute the matrix <& = E^Es-f, where E s -f and E s j. denote the submatrices extracted from E s by removing its first 
row and its last row, respectively. 

(v) Compute the eigenvalues of A;, i = 1,2,..., K T . 

(vi) Recover the unknown delays as follows: n = — ^arg(Ai). 



and recall that the ith element of a[n] is given by ©. We can therefore write the iV-lengfh sequence 
{di [n] } for each index i in the following matrix-vector form 

aj = XR(i/,)a„ (27) 

where is a length- N vector whose nth element is ai[n], X is an N x N diagonal matrix whose 
nth diagonal element is given by x n , R(fj) is an N x K v> i Vandermonde matrix with (n,j)th element 
e j2nu i:j nT^ an( j a . j g l en gQi-j( ui vector with jth element a-ij. The matrix X in (|27T ) can be inverted under 
the assumption that the sequence {x n } satisfies \x n \ > for every n = 0, . . . , N — 1, resulting in 

hi = R(v>i)ai, (28) 

where a^ = X _1 aj. From a simple inspection, we can express the elements of a^ as 

Oi[n] = " v '-' 2;a " " T . < n < N - 1. (29) 
i=i 

It is now easy to see from this representation that recovery of the Doppler-shifts from the sequences 
{oj[n]} is equivalent to the problem of recovering distinct frequencies from a (weighted) sum of complex 
exponentials. In the context of our problem, for each fixed index i, the frequency of the jth exponential 
is ujy = 2iTVijnT and its amplitude is aij. 

Fortunately, the problem of recovering frequencies from a sum of complex exponentials has been 
studied extensively in the literature and various strategies exist for solving this problem (see Ifl4l for a 
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review). One of these techniques that has gained interest recently, especially in the literature on finite rate 
of innovation |[30l - |[33l . is the annihilating-filter method. The annihilating-filter approach, in contrast to 
some of the other techniques, allows the recovery of the frequencies associated with the zth index even at 
the critical value of N = 2K V ^. On the other hand, subspace methods such as ESPRIT liT2l . matrix-pencil 
algorithm |fT3l , and the Tufts and Kumaresan approach ifTTIl are generally more robust to noise but also 
require more samples than 2K V ^. Once the Doppler-shifts for each index i have been recovered then, 
since RJ(fj)R(i/j) = I because of the requirement that N > 2K V ^, the attenuation factors {a^} are 
determined as 

oti = B)(vi)ai, i = l,...,K T . (30) 

IV. Sufficient Conditions for Identification 

Our focus in Section [III] was on developing a recovery algorithm for the identification of ULSs. 
We now turn to specify conditions that guarantee that the proposed procedure recovers the set of triplets 
{{Ti, v%ji Oiij)} that describe the parametric ULS %. We present these requirements in terms of equivalent 
conditions on the time-bandwidth product TW of the input signal x(t). This is a natural way to describe 
the performance of system identification schemes since TW roughly defines the number of temporal 
degrees of freedom available for estimating % flU. 

To begin with, recall that the recovery stage involves first determining the unknown delays r from the 
set of equations given by (fl4l ) (cf. Section IIII-BI ). Therefore, to ensure that our algorithm successfully 
returns the parameters of T~L, we first need to provide conditions that guarantee a unique solution to (fT4l ). 
To facilitate the forthcoming analysis, we let d [A] = {d [n] , n £ and b [A] = {b [n] , n £ Z} denote 
the set of all vectors d [n] and b [n], respectively. Using this notation, we can rewrite (fT4l) as 

d[A] =N(r)b[A] . (31) 

We now leverage the analysis carried out in iflOl to provide sufficient conditions for a unique solution 
to ( f3Tb ; see ifTOl for a formal proof. 

Proposition 1. If Or, b [A] / 0) solves ( f3Tb and if 

p > 2K T - dim {span (b [A])) (32) 

then (r,b[A]j is the unique solution of ( 1311 ). Here, span (h [A]) is used to denote the subspace of 
minimal dimensions that contains b [A]. 
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Proposition Q] suggests that a unique solution to (f3TT > — and, by extension, unique recovery of the 
set of delays r — is guaranteed through a proper selection of the parameter p. In particular, since 
dim (span (b [A])) is a positive number in general, we have from Proposition Q] that p > 2K T is a 
sufficient condition for unique recovery of r and b [A]. From Condition Q] in Section [Till we have seen 
that the parameter p effectively determines the minimum bandwidth W of the prototype pulse [cf. (|22|)1. 
Combining the condition p > 2K T and the one obtained earlier in d22l ) leads to the following sufficient 
condition on the bandwidth of the input signal for unique recovery of r and b [A] : 

W > -^A (33) 

The second step in the recovery stage involves recovering the Doppler-shifts and attenuation factors 
(cf. Section ITlI-BI ). We now investigate the conditions required for unique recovery of the Doppler-shifts. 
Recall that the vectors b[n] and a[n] are related to each other by the invertible frequency relation (fT2l ). 
Therefore, since the diagonal matrix D [e^ wT , r) is invertible and completely specified by r, the condition 
given in (l33l also guarantees unique recovery of the vectors a[n] from b[n]. Further, it can be easily 
verified that the matrix R(fj) in (|28T ) has the same parametric structure as that required by Proposition Q] 
We can therefore once again appeal to the results of Proposition Q] in providing conditions for unique 
recovery of the Doppler-shifts {v{\ from the vectors {a^} [cf. (l28l)1. To that end, we interchange p with 
N and K T with K v ^ in Proposition Q] and use the fact that dim (span (oj)) = 1 (since it is a nonzero 
vector). Therefore, by making use of Proposition [Q we conclude that a sufficient condition for unique 
recovery of vi from (|28T ) is 

N > 2K Vli . (34) 



Condition (1341 ) is intuitive in the sense that there are 2K Vj i unknowns in (I28t (K u> i unknown Doppler- 
shifts and K V j unknown attenuation factors) and therefore at least 2K V ^ equations are required to solve 
for these unknown parameters. Finally, since we need to ensure unique recovery of the Doppler-shifts 
and attenuation factors for each distinct delay n, we have the condition 

N > 2 max K vi (35) 

i 

which trivially ensures that d34l holds for every i = 1,...,K T . We summarize these results in the 
following theorem. 

Theorem 2 (Sufficient Conditions for System Identification). Suppose that % is a parametric ULS that 
is completely described by a total of K = Y^=i^^,i triplets (fi^ij^ij)- Then, irrespective of the 
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distribution of {(t%, within the delay-Doppler space, the recovery procedure specified in Section 17771 
with samples taken at {t = 2nn/W} uniquely identifies % from a single observation %(x(t)) as long 
as the system satisfies [A1]—[A3], the probing sequence {x n } remains bounded away from zero in the 
sense that \x n \ > for every n = 0, . . . , N — 1, and the time— bandwidth product of the (known) input 
signal x(t) satisfies the condition 

TW > 87rK T K U)max (36) 

where K Ut7nax = maxj K v i is the maximum number of Doppler-shifts associated with any one of the 
distinct delays. 

Proof: Recall from the previous discussion that the delays, Doppler-shifts, and attenuation factors 
associated with H can be uniquely recovered as long as N > 2K u>max , W > 47r jT T , and p > 2K T . Now 
take N = and note that under the assumption TW > %ixK T K v ^ max , we trivially have N > 2K u ^ max . 
Further, since T = NT and since sampling rate of 2-k/W implies p = 2ir/(WT), we also have that 
W = ^Kr. ^ p = 2K T , completing the proof. ■ 
Theorem [2] implicitly assumes that K T (or an upper bound on K T ) and K v ^ max (or an upper bound 
on K v ^ max ) are known at the transmitter side. We explore this point in further detail in Section |V] and 
numerically study the effects of "model-order mismatch" on the robustness of the proposed recovery 
procedure. It is also instructive (especially for comparison purposes with related work such as O, ifTTl ) 
to present a weaker version of Theorem [2] that only requires knowledge of the total number of delay- 
Doppler pairs K. 

Corollary 1 (Weaker Sufficient Conditions for System Identification). Suppose that the assumptions of 
Theorem\2\hold. Then the recovery procedure specified in Section\III\with samples taken at {t = 2mr/W} 
uniquely identifies % from a single observation T~L(x(t)) as long as the time— bandwidth product of the 
known input signal x(t) satisfies the condition TW > 2n(K + l) 2 . 

Proof: This corollary is a simple consequence of Theorem |2] and the fact that K T K V)Tnax < ^ K+ ^ . 
To prove the latter fact, note that for any fixed K and K T , we always have K V)max < K — (K T — 1). 
Indeed, if K UyTnax were greater than K — (K T — 1) then either Y^n=i > K or there exists an i such 
that K Vf i = 0, both of which are contradictions. Consequently, for any fixed K, we have that 

K T K umax < -Kl + (K + 1)K T (37) 

and since the maximum of —A' 2 + (K + l)K T occurs at K T = -^f-, we get K T K UtTnax < ^ K+ ^ . ■ 
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V. Discussion 

In Sections [III] and UVJ we proposed and analyzed a polynomial-time recovery procedure that ensures 
identification of parametric ULSs under certain conditions. In particular, one of the key contributions of 
the preceding analysis is that it parlays a key sub-Nyquist sampling result of fTOl into conditions on the 
time-bandwidth product, TW, of the input signal x(t) that guarantee identification of arbitrary linear 
systems as long as they are sufficiently underspread. Specifically, in the parlance of system identification, 
Corollary Q] states that the recovery procedure of Section [III] achieves infinite simally -fine resolution in the 
delay-Doppler space as long as the temporal degrees of freedom available to excite a ULS are on the 
order of Q(K 2 ). In addition, we carry out extensive numerical experiments in Section [VTT1 which confirm 
that — as long as the condition TW > 2n(K + l) 2 is satisfied — the ability of the proposed procedure to 
distinguish between (resolve) closely spaced delay-Doppler pairs is primarily a function of the signal-to- 
noise ratio (SNR) and its performance degrades gracefully in the presence of noise. In order to best put 
the significance of our results into perspective, it is instructive to compare them with corresponding results 
in recent literature. We then discuss an application of these results to super-resolution target detection 
using radar in Section |VT] 

There exists a large body of existing work — especially in the communications and radar literature — 
treating identification of parametric ULSs; see, e.g., (2], J9], lfT5l - |[T9"1 . One of the approaches that is 
commonly taken in many of these works, such as in J9], lTT5l - |[T8l . is to quantize the delay-Doppler 
space (r, v) by assuming that both t\ and v%j lie on a grid. The following theorem is representative of 
some of the known results in this case[] 

Theorem 3 ((9], 111710 . Suppose that % is a parametric ULS that is completely described by a total 
of K = ^^Z\Kvi triplets (r^, Uij, c%). Further, let the delays and the Doppler-shifts of the system 
be such that ti = r^W -1 and = tiff -1 for rj G Z + and £ij £ Z. Then % can be identified in 
polynomial-time from a single observation ~H(x(t)) as long as the system satisfies [A1]—[A3] and the 
time-bandwidth product of the input signal x(t) satisfies TW = Q(K 2 / log TW). 

There are two conclusions that can be immediately drawn from Theorem [3] First, both |9], ifTTl require 
about the same scaling of the temporal degrees of freedom as that required by Corollary[T} TW ~ Q(K 2 ). 
Second, the resolution of the recovery procedures proposed in (9], IfTTl is limited to W _1 in the delay 

3 It is worth mentioning here that a somewhat similar result was also obtained independently in 1341 in an abstract setting. 
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space and T _1 in the Doppler space because of the assumption that n = r^W^ 1 and Uij = ^•7" -1 u 
Similarly, in another related recent paper lTT8l . two recovery procedures are proposed that have been 
numerically shown to uniquely identify % as long as TW S> 1 and each (r^ify) corresponds to one 
of the points in the quantized delay-Doppler space with resolution proportional to W^ 1 and T^ 1 in the 
delay space and the Doppler space, respectively. Note that the assumption of a quantized delay-Doppler 
space can have unintended consequences in certain applications and we carry out a detailed discussion 
of this issue in the next section in the context of radar target detection. 

Finally, the work in |[T9l leverages some of the results in DOA estimation to propose a scheme for 
the identification of linear systems of the form (O without requiring that r, = rjW -1 and = i^l 1 . 
Nevertheless, our results differ from those in |[T9ll in three important respects. First, we explicitly state the 
relationship between the time-bandwidth product TW of the input signal x(t) and the number of delay- 
Doppler pairs K = X)i=i Kv.i that guarantees recovery of the system response by studying the sampling 
and recovery stages of our proposed recovery procedure. On the other hand, the method proposed in fT9ll 
assumes the sampling stage to be given and, as such, fails to make explicit the connection between the 
time-bandwidth product of x(t) and the number of delay-Doppler pairs. Second, the algorithms proposed 
in (I9\i have exponential complexity, since they require exhaustive searches in a if -dimensional space, 
which can be computationally prohibitive for large-enough values of K. Last, but not the least, recovery 
methods proposed in |[T9ll are guaranteed to work as long as there are no more than two delay-Doppler 
pairs having the same delay, maxj K V) i < 2, and the system output is observed by an M-element antenna 
array with M ^ K. In contrast, our recovery algorithm does not impose any restrictions on the distribution 
of {(Tj, i^ij)} within the delay-Doppler space and is guaranteed to work with a single observation of the 
system output. 

VI. Application: Super-Resolution Radar 

We have established in Section [TV] that the polynomial-time recovery procedure of Section [Til] achieves 
infmitesimally-fine resolution in the delay-Doppler space under mild assumptions on the temporal degrees 
of freedom of the input signal. This makes the proposed algorithm extremely useful for application areas in 
which the system performance depends critically on the ability to resolve closely spaced delay-Doppler 

4 Note that there is also a Bayesian variant of Theorem [3] in (9| that requires TW ~ SI (A") under the assumption that H has a 
uniform statistical prior over the quantized delay-Doppler space. A somewhat similar Bayesian variant of Corollary [JJ can also 
be obtained by trivially extending the results of this paper to the case when H is assumed to have a uniform statistical prior 
over the non-quantized delay-Doppler space. 
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Fig. 3. Quantized representation of nine targets (represented by *) in the delay-Doppler space with T max = 10 /is and 
Vmax — 10 kHz. The quantized delay-Doppler approximation of the targets corresponds to W = 1.2 MHz and T = 0.48 ms. 

pairs. In particular, our method can be used for super-resolution target detection using radar. This is 
because the noiseless, clutter-free received signal in the case of monostatic radars is exactly of the form 
(Q]) with each triplet (r^, u^, a^) corresponding to an echo of the radar waveform x(t) from a distinct 
target ojf] The fact that our recovery procedure allows to identify arbitrary parametric ULSs, therefore, 
enables us to distinguish between multiple targets even if their radial positions are quite close to each 
other and/or their radial velocities are similar — the so-called super-resolution detection of targets. 

On the other hand, note that apart from the fact that none of the methods based on the assumption of a 
quantized delay-Doppler space can ever carry out super-resolution target detection, a major drawback of 
the radar target detection approach in works such as (9l, |fl"8l is that targets in the real-world do not in fact 
correspond to points in the quantized delay-Doppler space, which causes leakage of their energies in the 
quantized space. In order to elaborate further on this point, define L = \WT max ] and M = \Tv max j2\ 
and note that (canonical) quantization corresponds to transforming the C = [0, T max ] x [— u max /2, u max /2] 
continuous delay-Doppler space into a Q = {0, . . . , L} x {— M/2, . . . , M/2} two-dimensional quantized 

grid, which in turn transforms the received signal %(x(t)) at the radar into ll35l . ll36l Chapter 4] 

L M 

H(x(t))nJ2 E a l imx{t-fi)<* 7l ** mt (38) 

1=0 m=-M 

where ct£ m = Ylf=i aije J7r ( m-7 ~ l/ «)sinc(m— TVij) sine TWrj) and the quantized delay-Doppler 

5 In the radar literature, the term "monostatic" refers to the common scenario of the radar transmitter and the radar receiver 
being collocated. 



December 8, 2010 



DRAFT 



20 




• 


* 




True Targets 
Estimated Targets 






• 








• 
• 

• 


















• 






• 










• 






0.1 0.2 0.3 0.4 0.5 0.6 0.7 
Delay (x T max ) 

(b) 


8 0.9 



Fig. 4. Comparison between the target-detection performance of matched-filtering and our proposed recovery procedure for the 
case of nine targets (represented by *) in the delay-Doppler space with T max = 10 /is, v ma x = 10 kHz, VV = 1.2 MHz, and 
T = 0.48 ms. The sequence {x n } corresponds to a random binary (±1) sequence with iV = 48, the pulse g(t) is designed to 
have a nearly-flat frequency response in the working spectral band T , and the pulse repetition interval is taken to be T = 10 /is. 
(a) Target detection by matched-filtering the received signal H(x(t)) with the input signal x(t). (b) Target detection using the 
proposed recovery procedure with p = 12. 



pairs (?£, v m ) £ Q. It is now easy to conclude from (l38l) that, unless the original targets (delay-Doppler 
pairs) happen to lie in Q, most of the attenuation factors {aem} will be nonzero because of the sine 
kernels — the so-called "leakage effect." This has catastrophic implications for target detection using radar 
since leakage makes it impossible to reliably identify the original set of delays and Doppler-shifts. This 
limitation of target-detection methods that are based on the assumption of a quantized delay-Doppler 
space is also depicted in Fig. [3] for the case of nine hypothetical targets. The figure illustrates that each 
of the nine non-quantized targets not only contributes energy to its own (fi, v m ) in Q but also leaks its 
energy to the nearby points in the quantized space. 

Owing to the fact that leakage can cause missed detections and false alarms, conventional radar literature 
in fact tends to focus only on recovery procedures that do no impose any structure on the distribution 
of {{Ti,Uij)} within the delay-Doppler space. The most commonly used approach in the radar signal 
processing literature corresponds to matched-filtering (MF) the received signal with the input signal x(t) 
in the delay-Doppler space J2|- The MF output x(t, v) takes the form 

X (t, v) = / H(x(t))x*(t - t) exp(-j2<irvt)dt = J2J2 a ^ A{ - T ~ Ti > v ~ u v) (39) 

where A{t,v) = J R x(t)x*(t — r) exp(— j2-Kvt)dt is termed the Woodward's ambiguity function of 
x(t). It can be easily deduced from (|39l ) that the resolution of the MF-based recovery procedure is tied 
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Fig. 5. Delay-Doppler representation of a parametric ULS H corresponding to K — 6 delay-Doppler pairs with T max = 10 /xs 
and v m ax = 10 kHz. 

to the support of the ambiguity function in the delay-Doppler space. Ideally, one would like to have 
A(t, v) = 5{t)8{v) for super-resolution detection of targets but two fundamental properties of ambiguity 
functions, namely, |.A(0,0)| 2 = f\x(t)\ 2 dt and JJ \A(t, v)\ 2 dTdv = f \x(t)\ 2 dt, dictate that no real- 
world signal x(t) can yield infmitesimally-fine resolution in this case either 0. In fact, the resolution 
of MF-based recovery techniques also tends to be on the order of and T _1 in the delay space and 
the Doppler space, respectively, which severely limits their ability to distinguish between two closely- 
spaced targets in the delay-Doppler space. This inability of MF-based methods to resolve closely-spaced 
delay-Doppler pairs is depicted in Fig. 01 This figure compares the target-detection performance of MF 
and the recovery procedure proposed in this paper for the case of nine closely-spaced targets. It is easy 
to see from Fig. HJa) that matched-filtering the received signal %{x{t)) with the input signal x(t) gives 
rise to peaks that are not centered at the true targets for a majority of the targets. On the other hand, 
Fig. ffib) illustrates that our recovery procedure correctly identifies the locations of all nine of the targets 
in the delay-Doppler space. 

VII. Numerical Experiments 

In this section, we explore various issues using numerical experiments that were not treated theoretically 
earlier in the paper. These include robustness of our method in the presence of noise and the effects 
of truncated digital filters, use of finite number of samples, choice of probing sequence {x n }, and 
model-order mismatch on the recovery performance. Throughout this section, the numerical experiments 
correspond to a parametric ULS H that is described by a total of K = 6 delay-Doppler pairs with 
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K T = 2 and K v< \ = K Vj 2 = 3. The locations of these pairs in the delay-Doppler space are given by 
Fig. [5j while the attenuation factors associated with each of the six delay-Doppler pairs are taken to have 
unit amplitudes and random phases. 

In order to identify %, we design the prototype pulse g(t) to have a constant frequency response over 
the working spectral band T = [— ^p, yjp] with p = 4 and T = 10 /is, that is, G(u>) 1 when u G T 
and G(ui) « when u ^ T . In other words, the input signal x(t) is chosen to have bandwidth W = 
In addition, unless otherwise noted, we use a probing sequence {x n } corresponding to a random binary 
(±1) sequence with N = 30, which leads to a time-bandwidth product of TW « 2407T. Note that the 
chosen time-bandwidth product here is more than the lower bound of Theorem [2] by a factor of 5 so as 
to increase the robustness to noise. Also, unless otherwise stated, all experiments in the following use an 
ideal (flat) LPF as the sampling filter (cf. Fig. [2]). We use the ESPRIT method described in Section [III] 
for recovery of the delays and the matrix-pencil method |[T3l for recovery of the corresponding Doppler- 
shifts. Given the rich history of these two subspace methods, there exist many standard techniques in 
the literature (see, e.g., ||37l , (38}) for providing them with reliable estimates of the model orders in 
the presence of noise. As such, we assume in the following that both these methods have access to 
correct values of K T and K u /s. Finally, the performance metrics that we use in this section are the 
(normalized) mean-squared error (MSE) of the estimated delays and Doppler-shifts (averaged over 100 
noise realizations), defined as 

e delay = „E K fi ~ T i)/ T ™ax) 2 > ( 40 ) 
i=l 

and 

1 2 3 

e Doppler = g Yl £ IP*' ~ "ij)/»rnax] 2 , (41) 
i=l j=l 

where fj and j>y denote the estimated delays and Doppler-shifts, respectively. 

1) Robustness to Noise: We first examine the robustness of our method when the received signal 
7i(x(t)) is corrupted by additive noise. The results of this experiment are shown in Fig. [6j which plots 
the MSE of the estimated delays and Doppler-shifts as a function of the SNR. It can be seen from the 
figure that the ability of the proposed procedure to resolve delay-Doppler pairs is primarily a function 
of the SNR and its performance degrades gracefully in the presence of noise. 

2) Effects of Truncated Digital-Correction Filter Banks: Recall from Section ITTT1 that our recovery 
method is composed of various digital-correction stages (see also Fig. [2]). The filters used in these stages, 
which include {<fo[ra]} and {^[n]}, have infinite impulse responses in general so that their practical 
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Fig. 6. Mean-squared error of the estimated delays and Doppler-shifts as a function of the signal-to-noise ratio. 
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Fig. 7. Mean-squared error of the estimated delays and Doppler-shifts as a function of the signal-to-noise ratio for various 
lengths of the impulse responses of the filters. 



implementation requires truncation of their impulse responses. The truncated lengths of these filters also 
determine the (computational) delay and the computational load of the proposed procedure. Fig. [7] plots 
the MSE of the estimated delays (Fig. |7Ja)) and Doppler-shifts (Fig. E^b)) as a function of the SNR for 
various lengths of the impulse responses of the filters. There are two important insights that can be drawn 
from the results of this experiment. First, for a fixed length of the impulse responses, there is always some 
SNR beyond which the estimation error caused by the truncation of the impulse responses becomes more 
dominant than the error caused by the additive noise (as evident by the error floors in Fig. |7]). Second, 
and perhaps most importantly, filters with 49 taps seem to provide good estimation accuracy up to an 
SNR of 60 dB, whereas filters with even just 35 taps yield good estimates at SNRs below 50 dB. 
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3) Effects of Finite Number of Samples: The sampling filter used at the front-end in Fig. [2] is 
bandlimited in nature and, therefore, has infinite support in the time domain. Consequently, our sampling 
method theoretically requires collecting an infinite number of samples at the back-end of this filter. The 
next numerical experiment examines the effect of collecting a finite number of samples on the estimation 
performance. The results are reported in Fig. [U which depicts the MSE of the estimated delays (Fig. Ufa)) 
and Doppler-shifts (Fig. Hfb)) as a function of SNR for different numbers of samples collected at the 
output of the sampling filter (corresponding to an ideal LPF). As in the case of truncation of digital- 
correction filter banks, there is always some SNR for every fixed number of samples beyond which the 
estimation error caused by the finite number of samples becomes more dominant than the error due to 
additive noise. Equally importantly, however, note that x(t) in these experiments corresponds to a train 
of ./V = 30 prototype pulses. Therefore, under the assumption of p = 4 samples per pulse period T, it is 
clear that we require at least N • p = 120 samples in total to represent just the input signal x(t). On the 
other hand, Fig. [8] shows that collecting 248 samples, which is roughly twice the minimum number of 
samples required, provides good (delay and Doppler) estimation accuracy for SNRs up to 70 dB. 

Finally, it is worth noting here that making use of an ideal LPF as the sampling filter requires collecting 
relatively more samples at the filter back-end due to the slowly-decaying nature of the sine kernel. 
Therefore, in order to reduce the number of samples required at the back-end of the sampling filter for 
reasonable estimation accuracy, we can instead make use of sampling filters whose (time-domain) kernels 
decay faster than the sine kernel. One such possible choice is a raised-cosine filter with roll-off factor 
equal to 1, whose frequency response is given by S(u) = ^ ( 1 + cos(-^w) ) when uj £ T and S(uj) = 
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Fig. 9. Frequency response of a raised-cosine filter with roll-off factor 1. 
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Fig. 10. Mean-squared error of the estimated delays and Doppler- shifts as a function of the signal-to-noise ratio for different 
numbers of samples collected at the output of a raised-cosine sampling filter with roll-off factor 1. 



when oj T . It is a well-known fact (and can be easily checked) that this filter decays faster in the time 
domain than the sine kernel. However, the main issue here is that raised-cosine filter does not satisfy 
Condition [2] in Section [TTll since its frequency response is not bounded away from zero at the ends of 
the spectral band F (see, e.g., Fig. [9]). 

However, we now show that this problem can be overcome by slightly increasing the sampling rate 
and the bandwidth requirement stated in Section JV] Specifically, note that Proposition Q] requires that 
the parameter p, which controls the minimal bandwidth of x(t) and the sampling rate of our proposed 
procedure, satisfies p > 4 under the current simulation setup (since K T = 2). We now instead choose 
p = 6 and argue that raised-cosine filter can be successfully used under this choice of p. To this end, 
recall from Section [III] that the function of the digital-correction filters ipi [n] and tpQ [n] is to invert the 
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frequency response of the sampling kernel corresponding to the frequency bands denoted by 1 and 6 in 
Fig. respectively (under the assumption that the pulse g(t) has a flat frequency response). In the case 
of a raised-cosine filter, however, we cannot compensate for the non-flat nature of these two bands since 
they are not bounded away from zero. Nevertheless, because of the fact that we are using p = 6, we can 
simply disregard channels 1 and 6 after the first digital-correction stage and work with the rest of the 
four channels (2-4) only. We make use of this insight to repeat the last numerical experiment using a 
raised-cosine filter and report the results in Fig. [TOj It is easy to see from Fig. [10] that, despite increasing 
p to 6, raised-cosine filter performs better than an ideal LPF using fewer samples. 

4) Effects of the Probing Sequence: Theorem [2] in Section [TV] stipulates that the choice of the probing 
sequence {x n } has no impact on the noiseless performance of the proposed recovery procedure as long 
as \x n \ > for every n = 0, . . . , N — 1. However, it is quite expected that {x n } will have an effect on 
the performance in the presence of noise and implementation issues related to truncated digital filters 
and use of finite number of samples. The next experiment examines this effect for four different choices 
of binary probing sequences of length N = 32 that periodically alternate between +1 and —1 every r 
entries. The results are reported in Fig. [TT] which depicts the MSE of the estimated delays (Fig. [TTf a)) 
and Doppler-shifts (Fig. fTTT b)) as a function of the SNR for probing sequences with r = 1,2,4, and 
32. We can draw two immediate conclusions from observing the results of this experiment. First, faster 
alternating probing sequences (in other words, sequences with higher frequency content) appear to provide 
better resilience against the truncation of digital filters and the use of finite number of samples. Second, 
the effect of the choice of probing sequence is less pronounced at low SNRs, since the error due to noise 
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Fig. 12. Effects of model-order mismatch on the performance of the proposed recovery procedure corresponding to H with 
K = 14 delay-Doppler pairs. 



at low SNRs dominates the errors caused by other implementation imperfections. 

5) Effects of Model-Order Mismatch: Our final numerical experiment studies the situation where the 
conditions of Theorem [2] do not exactly hold. To this end, we simulate identification of a parametric ULS 
with K T = 4 delays. For the first 3 delays we take K vi = 2, i = 1, 2, 3, whereas we choose K v ± = 8 
for the last delay. Finally, we take the prototype pulse g(t) as at the start of this section (with bandwidth 
W = fj£), but we use a probing sequence {x n } corresponding to a random (±1) sequence with N = 8. 
Clearly, this does not satisfy the conditions of Theorem |2] because of the large number of Doppler-shifts 
associated with the last delay {K v ^ = 8). 

The results of this numerical experiment are reported in Fig. |T2j It can be easily seen from the figure 
that, despite the fact that x{t) does not satisfy the conditions of Theorem [2j our algorithm successfully 
recovers the first three delays and the corresponding Doppler-shifts. In addition, the fourth delay is 
correctly recovered but (as expected) the Doppler-shifts associated with the last delay are not properly 
identified. Note that in addition to demonstrating the robustness of our procedure in the presence of 
model-order mismatch, this experiment also highlights the advantage of the sequential nature of our 
approach where we first recover the delays and then estimate the Doppler-shifts and attenuation factors 
associated with the recovered delays. The main advantage of this being that if the input signal does not 
satisfy N > 2K y ^ for some i then recovery fails only for the Doppler-shifts associated with the zth delay. 
Moreover, recovery of the ith delay itself does not suffer from the mismodeling and it will be recovered 
correctly as long as the bandwidth of x(t) is not too small. 



December 8, 2010 



DRAFT 



28 



VIII. Conclusion 

In this paper, we revisited the problem of identification of parametric underspread linear systems that 
are completely described by a finite set of delays and Doppler-shifts. We established that sufficiently- 
underspread parametric linear systems are identifiable as long as the time-bandwidth product of the 
input signal is proportional to the square of the total number of delay-Doppler pairs. In addition, we 
concretely specified the nature of the input signal and the structure of a corresponding polynomial- 
time recovery procedure that enable identification of parametric underspread linear systems. Extensive 
simulation results confirm that — as long as the time-bandwidth product of the input signal satisfies the 
requisite conditions — the proposed recovery procedure is quite robust to noise and other implementation 
issues. This makes our algorithm extremely useful for application areas in which the system performance 
depends critically on the ability to resolve closely spaced delay-Doppler pairs. In particular, our proposed 
identification method can be used for super-resolution target detection using radar. 
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